Monte Carlo study of degenerate groundstates and residual entropy 
in a frustrated honeycomb lattice Ising model 
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We study a classical fully-frustrated honeycomb lattice Ising model using Markov chain Monte 
Carlo methods and exact calculations . The Hamiltonian realizes a degenerate ground state manifold 
of equal-energy states, where each hexagonal plaquette of the lattice has one and only one unsatisfied 
bond, with an extensive residual entropy that grows as the number of spins A''. Traditional single- 
spin flip Monte Carlo methods fail to sample all possible spin configurations in this ground state 
efficiently, due to their separation by large energy barriers. We develop a non-local "chain-flip" 
algorithm that solves this problem, and demonstrate its effectiveness on the Ising Hamiltonian with 
and without perturbative interactions. The two perturbations considered are a slightly weakened 
bond, and an external magnetic fleld h. For some cases, the chain-flip move is necessary for the 
simulation to flnd an ordered ground state. In the case of the magnetic field, two magnetized ground 
states with non-extensive entropy are found, and two special values of h exist where the residual 
entropy again becomes extensive, scaling proportional to A'' In where (j) is the golden ratio. 



I. INTRODUCTION 

Ising frustration is a common ingredient in spin models 
designed to search for and study exotic physics. The 
prototypical example is the well-known triangular lattice 
antifcrromagnetic (AFM) Ising model, 
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which admits a ground state without long-range or- 
der (with power-law spin correlations), where an exten- 
sive number of degenerate (equal-energj^ configurations 
causes a residual (T = 0) entropy [l|, This classi- 
cal "manifold" of ground states is the fertile foundation 
from which one expects novel or exotic order to spring. 
For example, the triangular lattice AFM Ising Hamil- 
tonian, perturbed with a quantum transverse field, un- 
dergoes order-by-disorder to realize a long-range ordered 
quantum dimer state If instead the perturbation is 
a nearest-neighbor in-plane ferromagnetic quantum ex- 
change, the system reveals an exotic "supcrsolid" phase 
with coexisting diagonal and off-diagonal long-range or- 
der [1,11,11. Ultimately, one would like to construct a 
spin Hamiltonian that is a true T = quantum param- 
agnet, spin- liquid Q, or resonating valence-bond phase 
[a Q . This could open a new window on our understand- 
ing of the world of deconfinement, quantum number frac- 
tionalization, and topological order p^ . a role in which 
frustrated interactions are sure to play. 

In the case of purely classical systems, perturbations 
to frustrated Ising Hamiltonian are of utmost importance 
to actual material physics, as is well documented in the 
spin ices (llj - rare-earth titanates that realize to a very 
close approximation Ising models on the frustrated py- 
rochlore lattice. In addition to the Ising exchange, in 




FIG. I: (Color online) A fully-frustrated honeycomb lat- 
tice. Single lines (labelled a) represent antiferromagnetic 
bond interactions, and double (blue) lines (labelled b) rep- 
resent ferromagnetic interactions. The lattice illustrated has 
N = 2 X L X L sites, with L — 4, and periodic boundary 
conditions. Dots illustrate one possible ground-state config- 
uration of the Ising model, with black dots representing 
spin-up, and empty sites spin-down. A (red) ex denotes each 
unsatisfied bond. 



these materials the dipolar interaction strength is signifi- 
cantly large, and has been shown to be a critical ingredi- 
ent in the realization of the spin ice state [H, [II as well 
as the prediction for long-range order [T^ . 

A special class of classical Ising model is of particular 
theoretical interest due to the ability to map their ground 
states to hard-core dimer models, in which dimers live on 
the bonds of the respective dual lattice [1, . The sim- 
plest frustrated dimer model is the classical triangular 
lattice dimer model [l6j : extensions of this prototypical 
example are known to harbor interesting physical phe- 
nomena. The goundstate of this model is one where each 
site of the triangular lattice has one and only one dimer 
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emanating from it. Although no long-range order exists 
in this ground state, when constrained to a torus the 
model admits configurations which r nay be categorized 
into four distinct topological sectors |15| . Quantum ex- 
tensions of this model promote, among other things, a 
short-ranged resonating-valence bond phase with decon- 
fined fractional excitations (spinous) [8|. The quantum 
dimer model on the triangular lattice has also been used 
to motivate the design of topologically protected qbits 
[l7t . Clearly, such models are ideal playgrounds for the 
study of the exotic physics mentioned above. 

The ground state of the classical dimer model on the 
triangular lattice maps to the so-called fully - frustrated 
(FF) honeycomb lattice Ising model [3, I^i that 

model, antiferromagnetically interacting spins are placed 
on the sites of a honeycomb lattice (labelled a in Fig. [T|), 
with the exception that one bond per hexagon has an 
exchange of opposite sign, i.e. a ferromagnetically inter- 
acting nearest-neighbor pair (labelled b): 

= J^(-l)''<-^>-'SfS'|. (2) 

Here, 5t^ij).b is a delta function, with value on bonds 
a and unity on bonds b. In this paper, we use the con- 
vention = ±1/2. The ground state of this model is 
one where each hexagon of the honeycomb lattice has one 
and only one unsatisfied bond. If this unsatisfied bond 
is mapped to represent a dimer on the dual (triangular) 
lattice, one immediately sees that the spin configurations 
in the ground state of the FF honeycomb lattice Ising 
model can be matched to the ground state of the classi- 
cal dimer model on the triangular lattice (a close-packed 
model of dimers with hard cores), which has one and 
only one dimer connected to each site [l^. This ground 
state does not have a long-range order; rather, it is an ex- 
tensive manifold of equal-energy disordered states, \yhich 
produce a residual entropy of 5 ~ 0.214 per spin [l6l.[T8j. 

We note that the equivalency of the ground-states 
of the two models (FF honeycomb Ising and triangular 
dimer) is true for many choices of the pattern of the FM 
bonds (not just the one in Fig.[T]). In other words, there is 
a "gauge freedom" for which bonds we call ferromagnetic 
and which we call antiferromagnetic, as occurs with sim- 
ilar mixed-bond models Hi . Other patterns for the FM 
bonds could be chosen in Fig. [T]that have an equivalent 
ground state for the above Hamiltonian (or, the exchange 
of FM and AFM bonds, which is another gauge choice). 
Only when we explore perturbations to this Ising Hamil- 
tonian in Sections III CI and IIIII does our particular gauge 
choice become important, and we will discuss it more 
there. 

In this paper, we study the ground states of this model 
using an efficient Markov chain Monte Carlo (MCMC) 
method. As described in Section |TT1 conventional local 
(or "single-spin flip" (SSF)) MCMC algorithms fail to 
explore the entire degenerate ground state ergodically, 
which leads to incorrect simulation results, in particu- 
lar for perturbed models. The problem can be allevi- 



ated with global loop and cluster algorithms [19|, how- 
ever those designed for use on corner-sharing triangular 
or tetrahedral lattices [13, do not generalize to the FF 
honeycomb Ising model. Therefore, we develop a general 
chain-flip algorithm which allows for full ergodicity in 
the MCMC sampling of the ground-state manifold of the 
unperturbed model. We demonstrate how this restored 
ergodicity uncovers a phase transition to a long-range or- 
dered state in a perturbed model, whereas conventional 
algorithms with only local conflguration changes do not. 

In Section IIIII of this paper, we use our MCMC algo- 
rithm to explore the evolution of the ground state of the 
Hamiltonian in an applied external magnetic field, where 
we find two non-trivial higher-magnetization states with 
non-extensive entropy (scaling as \/N). In addition, for 
two critical field values bounding these states, we find 
special points of restored extensive entropy. Remarkably, 
one is able to calculate the values of these "reemergent" 
extensive entropies exactly, and we find that they scale 
proportion to N In 0, where <p is the golden ratio. These 
values are confirmed by our MCMC simulations. 



II. THE CHAIN-FLIP ALGORITHM 

The most powerful method to study classically frus- 
trated lattice spin models is MCMC. The design of effi- 
cient algorithms has made much of the past discoveries 
in the field possible, but it is not without its difficulties. 
In particular, the ability to study perturbed Hamiltoni- 
ans, and phenomena related to the lifting of macroscopic 
Ising ground-state degeneracies, is critically dependent 
on global "cluster" algorithms which are able to traverse 
degenerate (or nearly degenerate) configurations in order 
to discover energetically preferred ground states. In other 
words, local updates (like the SSF mentioned above) may 
suffer from a loss of ergodicity, or an exponential sup- 
pression of computational efficiency, effectively becoming 
frozen into particular states due to the presence of large 
energy barriers between configurations. However, global 
updates may restore ergodicity to the algorithm. This 
has been demonstrated with the development of loop al- 
gorithms in classical two-dimensional (2D) ice and vertex 
toy models [l^ , and has matured into a very general set 
of loop algorithms applicable to a wide range of models 
on corner-sharing triangular-based lattices in two dimen- 
sions (kagome) [20| and three dimensions (pyrochlore) 

A. Description of the Algorithm 

In this section, we examine in more detail the mech- 
anism by which SSF updates become inefficient in our 
FF honeycomb lattice Ising model Eq. ([2]), before devel- 
oping a global "chain-fiip" algorithm which restores er- 
godicity at low temperatures. Consider the classical SSF 
Metropolis algorithm, where configurational changes are 
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made by attempting to flip each spin individuaUy, find- 
ing the corresponding change in energy AE = iSafter — 
-^before; and then accepting the flip with probabihty 



In this as in all frustrated models, the SSF method works 
well at higher temperatures T, but as the temperature 
becomes lower, the system begins to "freeze" (or lose 
ergodicity) into its disordered manifold of equal-energy 
states. Each degenerate ground state configuration is at 
the bottom of a local energy well, which means that most 
"nearby" system configurations (configurations with only 
a few different spins) have a significantly higher energy. 
As a result, the probability P that any SSF which breaks 
out of the degenerate manifold is accepted becomes ex- 
ponentially low. To get from one ground state to another 
(i.e. to move from one energy minimum to an adjacent 
one), multiple consecutive single spin fiips are needed, 
and since the probability of any one spin being fiipped 
is low, the chance that multiple spins are flipped con- 
secutively is very small. This means that the simulation 
dynamics become effectively frozen, or non-ergodic. 

In order to overcome this difficulty, one requires an al- 
gorithm that flips multiple spins simultaneously in a way 
that bypasses the large energy barriers and tunnels from 
one ground state to the next. We achieve this by intro- 
ducing a "chain-flip" algorithm for the FF honeycomb 
Ising model. To understand how this chain move works, 
it is first useful to understand the structure of ground 
states in the model (see Fig. [T]for one example). In any 
given configuration, we call a bond unsatisfied if it has 
positive energy (-I-J/4), and satisfied otherwise (—J/4). 
The spin configurations that contribute to the ground 
state manifold are the ones in which every hexagon on 
the lattice has one and only one unsatisfied bond. If flip- 
ping a group of spins creates as many unsatisfied bonds 
as it removes, such a process is equivalent to moving from 
one degenerate ground state to another. This net zero- 
energy move is not possible with single-spin flips alone. 

Now we look at a method for finding a group of spins 
that creates as many unsatisfied bonds as it removes 
(Fig [2]) . We build up a chain of spins by selecting ver- 
tices one at a time and counting the number of unsatisfied 
bonds that fiipping the spin creates, versus how many it 
removes. Call the total number of unsatisfied bonds at 
any point in the algorithm the net unsatisfied bonds. It 
is useful to note that fiipping a spin changes the state 
of frustration of its three neighboring bonds. We note 
that if at any point in building the chain we encounter a 
hexagon not initially in its ground state (i.e. with more 
than one unsatisfied edge) we cancel the chain, thus en- 
suring the lattice remains locally in a ground state and 
ensuring detailed balance. 

Our first step is to pick a hexagon at random and la- 
bel its unsatisfied bond 61 (see Fig. [5]). Label the two 
hexagons that share that bond hi and h2. Then we ran- 
domly pick one of 61's spins, and label it si, and label 




FIG. 2: (Color online) Labeled Lattice 



the third hexagon that si is adjacent to /i3. Label the 
bond shared by /il and /i3 bond 52, and the bond shared 
by h2 and h3 bond 63. Note that since 62 and 63 are 
edges of hi and h2 respectively, neither of them can be 
unsatisfied. Now store si as the first spin in our chain, 
and fiip it. Flipping si makes 61 satisfied and 62 and 63 
unsatisfied, giving us -1-1 net unsatisfied bonds. 

We know that h3 initially had one unsatisfied bond, 
and /i3 has four bonds we have yet to consider: two ad- 
jacent to 62 and 63, and two opposite 62 and 63. Label 
the spins that si shares with 62 and 63 with s2 and s3 
respectively. If /i3's initial unsatisfied bond is adjacent 
to s2, then s2 is now adjacent to two unsatisfied bonds 
(its third bond, being part of hi, is satisfied). Thus we 
can flip s2 and create —1 unsatisfied bonds, giving our 
chain (si and s2) a total of net unsatisfied bonds when 
flipped, and we are done. This argument works similarly 
if ft,3's initial unsatisfled bond is adjacent to s3. So we 
see that the minimum number of spins that need to be 
flipped to complete the chain algorithm is two. Alterna- 
tively, the algorithm continues, if /i3's unsatisfled bond is 
opposite either 62 or 63. Incidentally, the basic two-spin 
chain flip is analogous to the elementary dimer plaquette 
moved generated by the kinetic term in typical quantum 
dimer models 

To explore the more general case, assume without loss 
of generality that the unsatisfled bond occurs opposite 
to 63. In this case, we want to flip s2, so label the sec- 
ond bond s2 shares with h3 bond 64, label s2's third 
bond 65, and label 64's other spin s4. Bonds 64 and 65 
are satisfied, so flipping s2 creates -f 1 unsatisfied bonds, 
for a total of +2 net unsatisfied bonds. Now label /i3's 
unsatisfied bond 66, and s4's third bond 67. Hexagon 
h3 shares 64 with hA, and label the hexagon adjacent to 
both h3 and hA hexagon hb. Spin s4 has two neighbor- 
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ing unsatisfied bonds, b6 and the newly unsatisfied &4, 
and one neigliboring satisfied bond, 67, wfiich we know 
is satisfied because it is part of h5, which already has 66 
as its unsatisfied bond. Thus flipping s4 creates —1 un- 
satisfied bonds, for a total of +1 net unsatisfied bonds. 
The chain must therefore continue; we now begin the re- 
cursive phase of the algorithm. 

Hexagon /i4 must initially have unsatisfied bonds, and 
so far we have looked at three of its bonds, 65, 64, and 67. 
Bond 65 and 67 have been made unsatisfied by spin flips 
we've already done. Label 65's second spin s5 and 67's 
second spin s6. If /i4's initial unsatisfied bond is adjacent 
to s5 we can flip s5, creating —1 unsatisfied bonds, for a 
total of net unsatisfied bonds. The algorithm is there- 
fore done. Similarly, if /i4's unsatisfied bond is adjacent 
to s6, one can fiip s6 and the chain is done. The recur- 
sive case comes when M's unsatisfied bond is opposite 64. 
Label the bond adjacent to s5 on hA bond 68, the bond 
adjacent to s6 on hi hexagon 69, and the bond opposite 
64 bond 610. We now randomly choose either 68 or 69 
and fiip both of its spins. Say without loss of generality 
that we chose 68. Since 65 and 610 arc unsatisfied and 
the other two bonds adjacent to 68 arc satisfied, flipping 
both of 68's spins results in unsatisfied edges, for a to- 
tal of +1 net unsatisfied bonds. Label the hexagon that 
shares 68 with h4 hexagon h6. 

The recursion comes by realizing that our current state 
is the same as it was before we flipped 68's two spins: h6 
has one unsatisfied edge that we have not found yet, on 
one of three sides that we have not considered. Further- 
more, two of /i6's spins have been flipped. Thus we can 
apply the same logic in choosing spins to add to the chain 
that we applied when choosing the last two. Specifically, 
if /i6's unsatisfied bond is the one opposite 68, then we 
continue recursing. If the unsatisfied bond is one of the 
other two possibilities, we can end the chain with net 
unsatisfied bonds, or in other words, with a chain that, 
when flipped, takes us to another ground state. 

In T = simulations of the Hamiltonian Eq. where 
the model is expected to be in the ground-state, spins on 
completed chains can be flipped with probability unity. 
However, in the case of perturbed Hamiltonians (see sec- 
tions HTC] and |TTT]), the energy change AE incurred by the 
proposed chain-flip must be calculated. Then, this pro- 
posed flip is accepted or rejected based on the Metropolis 
condition, Eq. ([3]). 



B. Chain moves in the unperturbed model 

In order to test the efficiency of the chain move in 
a real simulation, in this section we present results for 
two MCMC codes for the unperturbed model Eq. 
each employing a Metropolis algorithm with a Boltzmann 
probability Eq. ([3]). The first uses a conventional SSF 
Metropolis algorithm, the second a combination of SSF 
and chain-flip updates. Specifically, the MCMC step us- 
ing SSFs alone consists of attempting to flip each spin 



in the lattice twice. The MCMC step using chain moves 
consists of attempting SSF on each spin in the lattice 
once, followed by one attempted chain move for every 20 
spins in the lattice. This convention for MCMC steps 
was chosen so that the CPU time of the chain-flip as- 
sisted step is roughly equal to that of the step using SSFs 
alone, thereby allowing us to directly compare results 
without relying on formal autocorrelation measurements. 
We note of course that, in a working MCMC code, other 
conventions for the Monte Carlo "step" may be chosen 
by the practitioner. 

In this work, simulations were performed at finite tem- 
peratures and on finite system sizes ranging from a few 
hundred to tens of thousands of spins, using of order 10^ 
MCMC steps. From such simulations, we examine the 
impact of the chain move certain thermodynamic quan- 
tities, in particular, the energy E, specific heat C — 
dE/dT, and magnetization per spin M — ^/Nj^i^i- 
Further, we restrict ourselves to looking at two proce- 
dures for obtaining such finitc-T data: 1) an annealing 
(or slow-cooling) algorithm, and 2) a quenched (or rapid- 
cooling) algorithm, details of which are reported below. 

First; we examine data obtained through annealing 
simulations. An annealing procedure is often employed in 
simulations of models with long time-scales or glassy dy- 
namics, and is known to help reach even complex ground 
states using very simple local (SSF) algorithms. In an 
annealing run, the simulation is started at a high tem- 
perature T/J >> 1, and the usual MCMC algorithm 
(a series of equilibriation and production steps) is em- 
ployed. After sufficient data is gathered, the temperature 
is lowered by a small step, keeping the system configu- 
ration from the previous (higher-temperature) step. The 
MCMC algorithm is then repeated, and the temperature 
is lowered again until the system settles in its ground 
state. 

Figure [3] illustrates the results for the energy E, the 
specific heat C, and the magnetization M of the anneal- 
ing algorithm for a system of 20000 Ising spins in the 
Ising model. It is clear that the results for E and C are 
similar for both types of MCMC step, and that both re- 
alize the proper ground state of the model, with energy 
per spin E/N = — J/4. Further, integration of C/T (see 
e.g. Ref. [2l|) for both of these simulation runs reveals 
that the model retains a residual entropy in its ground 
state of S/N = 0.214 to within numerical (1%) accu- 
racy [ll|. The lack of difference between the SSF and 
the chain algorithm data can be explained largely as the 
success of annealing: even without the chain move, an- 
nealing allows the SSF algorithm alone to find a ground 
state. However, the single spin flips cannot move between 
degenerate ground states at very low T; E and C are sim- 
ply unaffected as every ground state has the same energy. 
This is not the case with the magnetization M, as seen 
in Figure [3l Clearly, the expectation that the ground 
state magnetization per spin should be tightly distributed 
around a mean of M = can be violated in the SSF 
algorithm, where degenerate configurations with higher 
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FIG. 3: (Color online) The energy, specific heat, and mag- 
netization of a simulation of 20000 Ising spins on the unper- 
turbed FF honeycomb model. 



magnetization can be frozen in, as illustrated. However, 
with the chain-flip algorithm, the expected convergence 
to M = is found with high accuracy. 

Our next observation is of the acceptance rate of the 
chain moves (Fig. [J). The single spin flips work well un- 
til lower temperatures are reached, at which point their 




FIG. 4: (Color online) The acceptance rate of single spin flips 
and chain moves across a range of temperatures. 
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FIG. 5: (Color online) The distribution of chain lengths on 
lattice of 20000 spins. 



acceptance rate drops off to zero. Comparing to Fig. [31 
this happens at a temperature where the system has re- 
alized the degenerate ground state. The chain moves, in 
contrast, have a very low acceptance rate at high temper- 
atures. This is due to the large number of chains being 
aborted during generation, when the construction algo- 
rithm encounters a hexagon with the "incorrect" number 
of unsatisfied bonds (i.e. not one). As the temperature 
lowers and the system reaches the ground state manifold, 
the chain-move acceptance rate climbs to unity. Clearly, 
a combination of SSF and chain-flips is needed in order 
to give a MCMC step with a reasonable acceptance rate 
across all temperatures. 

It is interesting to consider the size of chains produced 
by the algorithm. Note flrst that every chain must have 
an even size to ensure an even number of "boundary" 
bonds (see Fig. [2]). In Fig.[5l we see that in a histogram 
of data collected for several parameter values, the most 
common chain size is two, with the distribution of chains 
decreasing rapidly with their size. At high temperatures 
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chain-flips is virtually unaffected by the energy barriers, 
and finds the ground state in less than 50 MCMC steps. 
This figure is therefore a testament to the increased efh- 
ciency of the algorithm with the chain move in reaching 
the ground state. 

In this section, we have demonstrated that MCMC 
simulations employing the chain-flip algorithm can both 
realize the ground state of the unperturbed FF honey- 
comb lattice Ising model more efficiently that those em- 
ploying single spin flips alone, and also remain unfrozen 
once in this ground state. We now demonstrate the use- 
fulness of the chain move in simulations where the Hamil- 
tonian has been perturbed by a small interaction that 
lifts the degeneracy of the ground state manifold. 



FIG. 6: (Color online) The energy of a simulation of 20000 
Ising spin in the unperturbed model at T / J = 0.05, as a 
function of Monte Carlo (MC) step. 



this distribution is weighted heavily towards chains of 
size two, but as temperature decreases the distribution 
flattens out somewhat, making longer chains more prob- 
able. We also note that the distribution of chain sizes 
is essentially unaffected by lattice size, since apparently 
very few chains are long enough to act on more than a 
very local area of the lattice. 

We complete our examination of the two possible 
MCMC algorithms mentioned above, turning now to a 
discussion of the quenching algorithm. Quenching refers 
to the procedure whereby a simulation is run completely 
at a single (usually low) temperature, in an attempt to 
obtain the ground state quickly without annealing. The 
desire to use a quenching algorithm is obvious if one is 
interested only in the low-temperature properties of the 
model, which is often the case, and it is widely used since 
it is also the easier of the two algorithms to implement. 
However, without the history of higher-temperature con- 
figurations provided by an annealing procedure, it is of- 
ten observed that simulations at low temperature have a 
more difficult time settling into their true ground state 
configuration, as it is possible to get trapped in local 
energy minima separated by large energy barriers that 
single-spin flips have difficulty overcoming. 

Figure [5] illustrates the energy of three different 
quenched simulations, all begun in a random initial state 
(at step zero), as a function of the number of MCMC 
steps. It is clear that the chain-flip simulation reaches 
the proper ground state energy in a fraction of the num- 
ber of steps that the simulations without it take. The 
effects of the energy barriers and local minima on simu- 
lations using only single spin flips is most obvious in the 
plot of the energy of the smaller lattice of size N = 1152. 
The plateaus correspond to local minima in the energy, 
and the jumps between plateaus correspond to several 
spins being flipped consecutively, allowing the simulation 
to "bypass" the energy barrier. The simulation using 



C. Chain moves in a perturbed model 

In physical cases of interest, for example in the model- 
ing of real materials, one typically expects a more com- 
plicated Hamiltonian than Eq. often taking the form 
of small perturbations added to (or modifying) the sim- 
ple Ising interaction. In many applications, these per- 
turbations require no modification of existing MCMC 
schemes. However in the case where the unperturbed 
ground state is an extensively degenerate manifold of 
equal-energy states (such as is the case with the FF hon- 
eycomb Ising model), this is not true. Specifically, if 
small perturbative interactions lift the degeneracy of the 
manifold by energetically favoring one or more specific 
configurations (e.g. promoting long-range order), it has 
been demonstrated that SSFs in a MCMC scheme can 
fail to find this true ground state [2l| . 

In our FF honeycomb model, one may see the dynami- 
cal freezing of the SSFs by inspecting the acceptance rate 
in Fig. m which also suggests that the chain-flips success- 
fully explore the degenerate manifold of states at very low 
temperatures. We test this idea by introducing a small 
perturbation in the Hamiltonian which slightly weakens 
the ferromagnetic bonds in Fig. [TJ The Hamiltonian can 
be written as 

H = J^^S-Sj5(^ij')^a ^ J''^ S-SjSf^ij-j,, (4) 

(iJ) {ij} 

where the first term is the Hamiltonian for the antifer- 
romagnetic bonds, and the second for the ferromagnetic 
bonds (i.e. both J and J' are positive constants - see 
Fig. [T] for the bond labels a and b). In the following 
discussion, we set J' /J = 0.90. With this slight pertur- 
bation, we expect that the system will select a unique 
ground state with long-range order from the extensive 
manifold of states in the unperturbed model. Since the 
ferromagnetic bond is slightly weakened, this order will 
be one where the unsatisfied bond (one per hexagon) is 
placed uniquely on the ferromagnetic bond 6. We expect 
the energy per spin of this ground state to be — 41J/160, 
which is less than the — J/4 of the unperturbed model. 
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FIG. 7: (Color online) The energy and specific heat of a sim- 
ulation of 20000 Ising spins on the perturbed FF honeycomb 
model, Eq. (g)). 



Indeed Fig. [7] shows that the chain moves significantly 
modify the behavior of the MCMC simulation. Results 
presented there are for an annealing algorithm, and in 
contrast to the unperturbed case, it is clear that the 
MCMC using SSFs docs not find the correct ground state, 
while the MCMC using chain-flips does. In addition, the 
specific heat curve in Fig. [7] shows a large peak in the 
algorithm using the chain move that is not present above 
the noise seen in the algorithm not using the chain move. 
One observes a dynamical freezing of the spin configura- 
tion in the SSF algorithm, where the MCMC simulation 
no longer is able to sample low-lying states, and hence 
eventually freezes into a disordered state. In contrast, 
the chain algorithm is able to find a phase transition to 
a long-range ordered state, promoted by the perturbed 
Hamiltonian, as evident from the peak in C. Integration 
of this specific heat peak (over T) in Fig. [7] with the chain 
moves finds all of the expected ln(2) entropy, confirming 
the development of a unique groundstate at T = 0. In 
contrast, with the SSF algorithm only, the the full ln(2) 
entropy is not recovered by the integration, indicating 
that a long-range ordered ground state is not found. 



FIG. 8: (Color online) The magnetization per spin of a sim- 
ulation of 288 Ising spins {L = 12) in the model Eq. ((S} at 
T/J = 0.005. The plateaus at M = 1/8 and M = 1/4 corre- 
spond to partially-ordered states, described in the text. The 
(asymptotic) values of the residual T = entropy S of the 
ground state configurations are also labelled. The extensive 
entropy spikes which occur a the special values of h/J = 1 
and h/J = 3/2 are discussed in Section UlIB I 



III. GROUNDSTATES IN AN EXTERNAL 
MAGNETIC FIELD 

We now turn to a consideration of the fully-frustrated 
honeycomb lattice Ising model, with the physically im- 
portant perturbation of an external magnetic field: 



H 



4)*<'^>'''S'fS'| + /i^S'f 



(5) 



The inclusion of a symmetry-breaking magnetic field is 
known to lift the degeneracy of the ground-state mani- 
fold in some models, e.g. the frustrated Ising AFM on a 
triangular lattice [22l |. However, in other cases, such as 
the frustrated Ising AFM on the kagome lattice, this is 
not the case, and the presence of a perturbative external 
magnetic field reduces the degeneracy but does not lift it 
all together 



A. Magnetization plateaus and entropy spikes 

Using MCMC simulations with the SSF and chain-flip 
algorithm, wc study the ground state of the Hamiltonian, 
Eq. on the honeycomb lattice. Figure [S] shows the 
evolution of the magnetization per spin M as a function 
of the applied field. Immediately upon application, the 
field promotes the development of a plateau with mag- 
netization M = 1/8. Inspection of simulation configu- 
rations on this plateau reveal that it corresponds to a 
partial ordering of spins. 

This order is illustrated in Fig. [SJ it can be described 
in terms of horizontal zig-zag "rows" (labelled a to c? 
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in Fig. [9]). There, every second row with ferromag- 
netic bonds (labelled a) is fully polarized (all spin- up). 
This choice of ordering pattern is obviously not gauge- 
invariant, since the pattern of FM bonds has been cho- 
sen to break the lattice rotational symmetry. It can be 
seen however that it lowers the energy of the spins as- 
sociated with these bonds in the magnetic field, while 
retaining a configuration that is a member of the ground 
state manifold in the unperturbed model (i.e. one and 
only one unsatisfied bond per hexagon). The rows adja- 
cent to the fully polarized row, (i.e. every row without 
a ferromagnetic bond, or rows b and d in Fig [S]), are 
forced into a specific configuration, with alternating up 
and down spins, in order to maintain one unsatisfied bond 
per hexagon. Finally, the remaining rows with ferromag- 
netic bonds (the c rows in Fig [5]) each have two possible 
configurations, each alternating two up spins, with two 
down spins along their ferromagnetic bonds. Remark- 
ably, this is not a unique long-range ordered state in two 
dimensions, since the FM rows of spins c sandwiched in 
between the fully polarized FM rows a are ordered in- 
dependent of the FM c rows above and below. Thus, a 
ground state entropy remains. However, it is no longer 
extensive (i.e. scaling as L^), but scales as L, the lattice 
linear dimension. See Section IIIIBI below for an exact 
expression for this entropy. 




FIG. 9: (Color online) The ground-state spin configuration 
of the Hamiltonian Eq. ((5} at magnetization M = 1/8 (0 < 
h/J < 1). Black dots represent spin-up, and empty sites 
spin-down. The 16 sites making up the unit cell are outlined 
by the dashed (red) line. To obtain the M = 1/4 state, the 
down-spins on the zigzag "row" labelled c must be flipped to 
up-spins. 

The M = 1/8 semi-ordered state remains stable to 
moderately large applied fields (as evident by the plateau 
in Fig. [51), until for h/J> 1 a second plateau is reached 
at M = 1/4. We note that for this plateau, the ground 
state is forced out of the degenerate manifold of states 
(with one and only one unsatisfied bond per hexagon). 
Hence, chain moves cease to be effective (although SSFs 
still contribute), resulting in the increased noise in the 
magnetization plateau in Fig. [8] and a rounding of the 
associated transition. From observation of simulation 
configurations, it is apparent that this plateau corre- 
sponds to the flipping of all down-spins to up-spins on 



row c (Fig. [9]). This configuration is again not truly long- 
ranged ordered in two dimensions. Rather, in this case 
the antiferromagnetic rows b and d have become inde- 
pendent of each other, and the ferromagnetic row c has 
become fixed. Thus, we find that the entropy of the 1/4 
plateau is actually twice as large as the entropy of the 
1/8 plateau, although still scaling as the linear system 
size L (see Section IHI Bp . Finally, for applied magnetic 
field h/J > 3/2, the remaining down-spins on the b and 
d rows flip to up-spins, and the system becomes fully- 
polarized. 

A curious phenomenon occurs in the transition regions 
between the three magnetization plateaus illustrated in 
Fig. [51 For the special values h/J = 1 and h/J = 3/2, 
the model transitions from a configuration of decoupled 
quasi- ID ordered chains to once again being a disordered 
2D system. In other words, at precisely these critical 
values oi h/J, the model becomes "accidentally" macro- 
scopically degenerate due to fine-tuning of the magnetic 
field. At h/J = 1, this degeneracy occurs between states 
like in Fig. [9l and states where all spins on all zig-zag 
row c are up. Similarly, the degeneracy at h/J = 3/2 
occurs between states like this last state, and the M ~ 1 
state where the remaining spins (on the zig-zag rows b 
and d) flip up. For these two special field values where 
the fine-tuning of h/J causes an accidental degeneracy, 
we expect a reemergence of a residual ground state en- 
tropy, scaling as the system size N (similar to the h = 
case). Remarkably, one is able to calculate the values 
of these reemergent extensive entropies exactly in this 
model. This is shown in the next section, where we dis- 
cover that the asymptotic residual entropies are 

S'/A^ = In </-/8w 0.06015 for h/J^l, (6) 
S'/A^ = In 0/2 « 0.24061 for h/J ^3/2. (7) 

Here, « 1.618 is the golden ratio. Numerically, one 
can measure the values of the residual entropy via our 
MCMC simulations as the difference of the integral of 
C/T (over all T) from ln(2). For a lattice of size L = 32, 
we obtain 5* = 0.0604(2) at h/J = 1, and 5* = 0.2398(2) 
at h/J = 3/2, where finite-size trends clearly suggest 
that the MCMC approaches the above asymptotic results 
in the thermodynamic limit, to within error bars. The 
derivation of the exact asymptotic results is presented in 
detail in the next section. 



B. Exact entropy calculations 

In this section we derive exact expressions for the 
finite-size entropy of the system with external magnetic 
field for varying h/J. In our derivation, we concentrate 
on linear lattice sizes L that are a multiple of 4, based 
on the assumption (e.g. from Fig. [9] and discussions in 
Ref. 0) that the smallest unit cell that contains the 
ordered or partially-ordered structure is commensurate 
with L = 4. We also give approximate expressions for 
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large lattice size L, and asymptotic values for L tending 
to infinity. In our presentation, we refer to Fig. [5) 

1. <h/J <1 

As discussed above, for < h/J< 1, ground states are 
such that rows c can each have two configurations (with 
pairs of spins alternating up and down), and there are 
L/4 rows c. Rows a have all spins up, and rows b and d 
are fixed as in Fig. [51 This gives 2^/^ configurations, but 
the roles of rows a and c can be switched, which leads to 
an additional doubling of this number of configurations. 
In this way we obtain 

O = 2^/4+1, (8) 
5 = lnf7= (^^ + I^ln2, (9) 

— =(— + ^]\n2, (10) 
giving 5" ln(2)/8 • L w 0.087i in the limit of large L. 




FIG. 10: (Color online) The ground state spin configuration 
for h/J = 1. The binary digits 1 or are associated with rows 
of type B. Pairs of spins on rows of type A can be flipped if 
that row is "above" a B row of type 1, and "below" a B row 
of type 0, for example the spins circled by the dashed (red) 
line in row A2. Spins on rows Ai and A3 may not be flipped. 



4. h/J^l 



2. Kh/J < 3/2 

We now extend our notation and call rows a and rows 
c rows of type A (A = {a, c}) and, similarly, B = {6, d}. 
As discussed above, for 1 < h/J < 3/2, ground states 
arc such that each row B has its spins alternating up 
and down, and is thus in one of two possible states, inde- 
pendently from the other rows B: either all its up spins 
are adjacent to the row A above it, or to the row A below 
it. For each row B, let's assign a binary digit 1 to the 
former case (spins up are adjacent to the row A above 
the row B), and in the other case. Let m = L/2 be 
the number of rows B on the lattice. Then there are 2™ 
configurations for the rows _B, and wc get 



O = 2^/2 



L 



S = Infi = -ln2, 
2 



■In 2, 



S__ 

i.e. S w 0.173L in the limit of large L. 



(11) 
(12) 

(13) 



3. h/J > 3/2 

For h/J > 3/2 there is a single ground state, with all 
spins up, and we get 



= 1, 

S = Inn = 0, 

s 



N 



0. 



(14) 
(15) 

(16) 



We first observe that all ground states for 1 < h/J < 
3/2 are still ground states for h/ J = 1. (Recall that these 
states have all spins up on rows A, and spins alternating 
up and down on rows B). However, there are now many 
additional ground states, because, for every of the 2™ 
configurations of rows B, some of the rows A are allowed 
to flip some of their pairs of spins that are connected by 
a ferromagnetic bond from up to down. Indeed, rows A 
that have their upper neighbor row B in the state, and 
their lower neighbor row B in the 1 state, are allowed to 
flip pairs of spins without changing the energy, as long as 
no two adjacent pairs have spin down (see Fig. llOp . Note 
that this corresponds to the case that, for a pair of up 
spins in rows A, all neighboring spins are up as well. For 
these rows A, flipping a spin pair (which has neighboring 
pairs up) from up to down results in a magnetic energy 
change of 4 /i/2 (two spins flipped from up to down), and 
the energy change from the bonds is —4 J/2, because 
four bonds become satisfied. At h/ J — 1 the energy thus 
remains unchanged. We now want to count how many 
states can be obtained in this way. 

To this end, we first investigate how many rows A are 
allowed to flip pairs of spins, for a given configuration 
of rows B. Every configuration of rows B can be rep- 
resented by an m-digit binary number (see above), and 
every transition from a to a 1 in this binary number 
corresponds to a row A whose spin pairs can be flipped. 
(Note that wc have to include the periodic case, where 
the last digit is and the flrst digit is 1, in our count.) 
Let i?(m, J) be the number of m-digit binary numbers 
with d transitions from to 1. It is shown in Appendix 
A that 



i?(m,d) = 2 



m 
2d 



(17) 
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Next we have to count, for each row A that is allowed 
to flip its pairs, how many configurations there are in 
which no two adjacent pairs are down. Let n ~ L he 
the number of pairs on a row A, and let g{n) be the 
number of configurations in which no two adjacent pairs 
are down. (Again, we have to include the periodic case 
in our count.) Appendix B shows that g{n) satisfies the 
Fibonacci recurrence equation, with exact solution 



where 



Cl 



C2 



1.6180, 
-0.6180, 
V 1.1708, 
=i -0.1708. 



1 


+ V5 




2 


1 


-x/5 




2 


5 


+ 3V5 




10 


5 


-3\/5 


10 



(18) 

(19) 
(20) 
(21) 
(22) 



Note that, for large L, g{L) can be approximated well by 
g(L) = ci(/.^. (23) 
The number of ground states, O, is now given by 



L/4 

n = J2RiL/2,d)g{Lf. 



(24) 



d=0 



(Note that there arc at most d ^ L/4 transitions from 
to 1 in a binary number with m — L/2 digits.) 

It is shown in Appendix C that a closed-form expres- 
sion for this sum is given by 



L/2 



1 + V9{L)) +[l-V9{L) 



L/2 



and then 5 = In 17 and S/N = \nn/(2L^). 
For large L, we can approximate Eq. (pS)) as 



n^2g{L) 



L/4 



2ci 



L/4 ,lV4 



which leads to 



S 

S_ 

N 



Inn : 
In 2 



In 2 
1 

2L2 + ^ 



L 
In Cl 



In Cl 
1 

+ 8 



I? 
4 



ln( 



Inc 



(25) 

(26) 

(27) 
(28) 



We thus find an asymptotic value for S/N equal to 
In (A/8 « 0.06015. 



5. h/J = 3/2 

We saw above that, for 1 < h/J < 3/2, ground states 
have all spins up in all rows A, and spins alternating up 




FIG. 11: (Color online) The ground state spin configuration 
for h/J — 3/2. Spins on rows of type B can be flipped as 
long as the condition of not having two adjacent down spins 
is satisfied. Down spins circled by a dashed (red) line can be 
flipped to up spins without a change in energy. 



and down in rows B. These states are all ground states 
at h/J = 3/2 as well, but there are many additional 
ground states, because rows B can now flip some of their 
spins up or down without changing the energy, as long 
as no two adjacent spins are down (see Fig.[TT]). Indeed, 
flipping a spin (which has its three neighboring spins up) 
from up to down results in a magnetic energy change of 
2h/2 (one spin flipped from up to down), and the energy 
change from the bonds is —3 J/2, because three bonds 
become satisfied. At h/J = 3/2 the energy thus remains 
unchanged. We now want to count how many states can 
be obtained in this way. 

First, there are m = L/2 rows B that can change some 
of their spins independently. Second, every row B has 
n = 2 L spins that can be flipped up our down as long 
as no two adjacent spins are down. The number of valid 
configurations for each row B is thus given by g(2L), 
with g{n) the specific solution of the Fibonacci equation 
given in Eq. p8p . This gives 



n^g{nr^g{2Lf/^ 

and then S = Infi and S/N = \nn/{2L^). 
For large L, we can approximate Eq. (p9)) as 



n = g{2L) 



which leads to 



S- 

§_ 

N ■ 



Inn 
1 

II 



_ L 
~ 2 

In Cl 



Inci 
1 



ln< 



ln( 



(29) 

(30) 

(31) 
(32) 



We thus find an asymptotic value for S/N equal to 
\n<j)/2 w 0.24061. 



IV. DISCUSSION 

In this paper, we have developed a global chain-flip 
algorithm for Markov Chain Monte Carlo simulations 
of the fully frustrated honeycomb lattice Ising model. 
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Chain-flips are used to complement eonventional single- 
spin flips, in parameter regimes where the MCMC sim- 
ulation is required to explore the model's extensively- 
degenerate ground state manifold of minimally-frustrated 
spin configurations. We have demonstrated, through 
careful numerical simulations, that the chain-flip algo- 
rithm both increases simulation efficiency, and restores 
ergodicity in the sampling of the degenerate manifold of 
states. We have emphasized this latter point by demon- 
strating that chain-flips are necessary for the MCMC sim- 
ulation to find the proper ground state in the case where 
one of the members of the extensive manifold is made 
to have lower energy. In this perturbed model, chain- 
flips promote a low temperature phase transition to a 
long-range ordered state, that recovers all of the residual 
entropy of the unperturbed model. 

We have also used our MCMC algorithm to study 
the physically important extension of the FF honeycomb 
Ising model, where an external magnetic field h is ap- 
plied. In this case, moderate values of h promote the 
realization of partially-ordered states, corresponding to 
magnetization plateaus with values of M = 1/8 and 
M = 1/4. The precise nature of the partially-ordered 
states is dependent on the geometry with which frus- 
tration is introduced into the original (unperturbed) FF 
honeycomb Ising model, and is not gauge-invariant in 
this sense. An interesting phenomenon that occurs is 
the reemergence of extensive entropy "spikes" at h val- 
ues bounding the M = 1/4 plateau. Using a proof based 
on the reduction of the configurational disorder down to 
a Fibonacci recurrence, we are able to show that the en- 
tropy of the two reemergent spikes is equal to 



S/N 
S/N 



In 
In 



0.06015 
0.24061 



for 
for 



h/J 
h/J 



1, 

3/2. 



where cj) w 1.618 is the golden ratio. MCMC simulations 
confirm these results to a high degree of accuracy. It is 
interesting to note that, in one case (for h/J = 3/2), 
the entropy spike S/N = 0.241 is actually greater than 
the residual entropy for /i = {S/N ^ 0.214). The 
phenomenon of a magnetization plateau bounded by ex- 
tensive entropy spikes has previously been seen to occur 
on several models in one p3| and two dimensions [2^ . 
and also in 3D spin ice systems in an applied field along 
the [111] crystallographic direction [20j . 
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APPENDIX A: NUMBER OF m-DIGIT BINARY 
NUMBERS WITH d TRANSITIONS FROM TO 

1. 

Let i?(m, d) be the number of (periodic) m-digit binary 
numbers with d transitions from to 1 . We want to show 
that 



i?(TO, d) 



TO 

2d 



(Al) 



In a periodic binary number of size m, there are to pos- 
sible places to switch digits. Choose 2d of these to places 
in order to get d transitions from to 1. Once the 2c? 
locations where digits switch arc chosen, the number can 
be formed in two different ways (zeros and ones can be 
switched). This leads directly to expression (jAip . 



APPENDIX B: NUMBER OF CONFIGURATIONS 

dip). 

Consider a (periodic) line with n spins. We want to 
show that g(n\ the number of configurations in which 
no two adjacent spins are down, satisfies the Fibonacci 
recurrence equation. 



g{n) = g{n - 1) +3(71- 2), 
with exact solution 



5(n) = ci 0"-f 



(Bl) 



(B2) 



where 



Cl 



C2 



l + \/5 ^ 
2 

1 - V5 ^ 
2 

5 + 3\/5 

10 
5-3^5 

10 



1.6180, 
-0.6180, 
=rf 1.1708, 
=rf -0.1708. 



(B3) 



Note that this also covers the case of a (periodic) line 
with n fixed spin pairs, in which no two adjacent spin 
pairs are down. 

To count the number of valid states, we introduce two 
sets. First, C/„ is the set of valid states on rows of size 
n such that the first spin in the row is an up spin. We 
would like to find |J7„| recursively. We divide f7„ into 
two categories, rows starting with two up spins and rows 
starting with an up spin and then a down spin. An up 
spin followed by any member of Vn-x is in the first cat- 
egory, and all members of the first category must be of 
that form. If the row starts with an up spin and then a 
down spin, then the third spin must be an up spin, so the 
second category is every row of the form up spin, down 
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spin, then a member of Un-2- So 

Un = (T +Un-l) U (t + i +C/n-2) 



^ |i7„| = |[/„_i| + |[/„_2|. (B4) 

Now our second set is r„, the set of all valid states 
on rows of size n. Clearly [/„ C T„, so we need only 
worry about rows that begin with a down spin. If a row 
begins with a down spin, its second spin must be an up 
spin. Furthermore, since the row loops around due to 
the periodicity of the lattice, its last spin must be an up 
spin. Thus rows that begin in a down spin, end in an up 
spin, and have a member of Un-2 in between to account 
for all remaining members of r„. So 

Tn^Un U (i +Un-2+ T) 

^ |T„| = |t/„| + |[/„_2| 

= |[/„_l| + \Un-2\ + |t/„-3| + Wn-A. by (jB4| 

= |r„_i| + |r„_2| 

With |r„| ~ g{n), we obtain the Fibonacci equation, 
Eq. (|Bip . for g{n). The general solution of this recur- 
rence is given by Eq. (jB2p , with 4> and 9 the roots of the 
characteristic polynomial of the equation. The base cases 
(by inspection) are g{l) ~ 2 and g{2) = 3, which provide 
the values for the constants ci and C2 that are given in 
(|B3|) . 

APPENDIX C: CLOSED-FORM EXPRESSION 
FOR EQ. (HH). 

We show that 

m/2 

n = ^R{m,d)a'^, (CI) 

d=0 
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